The problem with GPs

Unless you are lucky (or clever), Gaussian process models are difficult to scale to large problems. For a Gaussian process \(\underset{n \times 1}{\boldsymbol{y}} \sim \mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})\):

A simple guide to computational complexity

How bad is the problem?

Practice - Migratory Model Prediction

After fitting the GP need to sample from the posterior predictive distribution at \(\sim3000\) locations \[ \boldsymbol{y}_{p} \sim \mathcal{N}\left(\mu_p + \Sigma_{po} \Sigma_o^{-1}(y_o - \mu_o) ,~ \Sigma_p - \Sigma_{po} \Sigma_{o}^{-1} \Sigma_{op}\right) \]

A bigger hammer?

Cholesky CPU vs GPU (P100)

Relative Performance

Aside (1) - Matrix Multiplication (P100)

Aside (2) - Memory Limitations

A general covariance is a dense \(n \times n\) matrix, meaning it will require \(n^2 \times\) 64-bits to store.

Other big hammers

More scalable solutions?

• Spectral domain / basis functions

• Covariance tapering

• GMRF approximations

• Low-rank approximations

• Nearest-neighbor models

Low rank approximations in general

Lets look at the example of the singular value decomposition of a matrix,

\[ \underset{n \times m}{M} = \underset{n \times n}{U}\,\underset{n \times m}{\text{diag}(S)}\,\underset{m \times m}{V^{\,t}} \]

where \(U\) are the left singular vectors, \(V\) the right singular vectors, and \(S\) the singular values. Usually the singular values and vectors are ordered such that the singular values are in descending order.

Example

Rank 2 approximation

Approximation Error

We can measure the error of the approximation using the Frobenius norm, \[ \lVert M-\tilde M\rVert_F = \left( \sum_{i=1}^m\sum_{j=1}^n (M_{ij}-\tilde M_{ij})^2\right)^{1/2} \]

Strong dependence

For a \(50 \times 50\) covariance matrix with a large effective range,

Weak dependence

For a \(50 \times 50\) covariance matrix with a short effective range,

How does this help?

There is an immensely useful linear algebra identity, the Sherman-Morrison-Woodbury formula, for the inverse (and determinant) of a decomposed matrix,

Aside - Determinant

Remember for any MVN distribution when evaluating the likelihood

we need the inverse of \(\Sigma\) as well as its determinant.

Low rank approximations for GPs

For a standard spatial random effects model,

\[ y(\boldsymbol{s}) = x(\boldsymbol{s}) \, \boldsymbol{\beta} + w(\boldsymbol{s}) + \epsilon, \quad \epsilon \sim N(0,~\tau^2 I) \]

\[ w(\boldsymbol{s}) \sim \mathcal{N}(0,~\boldsymbol{\Sigma}(\boldsymbol{s})), \quad \boldsymbol{\Sigma}(\boldsymbol{s},\boldsymbol{s}')=\sigma^2\;\rho(\boldsymbol{s},\boldsymbol{s}'|\theta) \]

if we can replace \(\boldsymbol{\Sigma}(\boldsymbol{s})\) with a low rank approximation of the form \(\boldsymbol{\Sigma}(\boldsymbol{s}) \approx \boldsymbol{U}\,\boldsymbol{S}\,\boldsymbol{U}^t\) where

• \(\boldsymbol{U}\) is \(n \times k\),

• \(\boldsymbol{S}\) is \(k \times k\), and

• \(A = \tau^2 I\) or a similar diagonal matrix

Gaussian Predictive Processes

For a rank \(k\) approximation,

• Pick \(k\) knot locations \(\boldsymbol{s}^\star\)

• PPs systematically underestimates variance (\(\sigma^2\)) and inflate \(\tau^2\), Modified predictive processs corrects this using

\[ \begin{aligned} \boldsymbol{\Sigma}(\boldsymbol{s}) \approx & \boldsymbol{\Sigma}(\boldsymbol{s},\boldsymbol{s}^\star) \, \boldsymbol{\Sigma}(\boldsymbol{s}^\star)^{-1} \, \boldsymbol{\Sigma}(\boldsymbol{s}^\star,\boldsymbol{s}) \\ &+ \text{diag}\Big(\boldsymbol{\Sigma}(\boldsymbol{s}) - \boldsymbol{\Sigma}(\boldsymbol{s},\boldsymbol{s}^\star) \, \boldsymbol{\Sigma}(\boldsymbol{s}^\star)^{-1} \, \boldsymbol{\Sigma}(\boldsymbol{s}^\star,\boldsymbol{s})\Big). \end{aligned} \]

Example

Below we have a surface generate from a squared exponential Gaussian Process where

Predictive Process Model Results

Performance

Parameter Estimates

Low Rank via Random Projections

Resulting approximation has a bounded expected error,

\[ E| \boldsymbol{A} - \boldsymbol{U}\boldsymbol{S}\boldsymbol{V}'\|_F \leq \left[1 + \frac{4\sqrt{k+p}}{p-1} \sqrt{\min(m,n)} \right] \sigma_{k+1}. \]

Random Matrix Low Rank Approxs and GPs

Low Rank Approximations and GPUs

Both predictive process and random matrix low rank approximations are good candidates for acceleration using GPUs.

• Both use Sherman-Woodbury-Morrison to calculate the inverse (involves matrix multiplication, addition, and a small matrix inverse).

• Predictive processes involves several covariance matrix calculations (knots and cross-covariance) and a small matrix inverse.

• Random matrix low rank approximations involves a large matrix multiplication (\(\boldsymbol{A}\,\boldsymbol{\Omega}\)) and several small matrix decompositions (QR, eigen).

Comparison \(n=15,000,\;k=\{100,\ldots,4900\}\)

Rand. Projection LR Depositions for Prediction

This approach can also be used for prediction, if we want to sample

\[ \boldsymbol{y} \sim \mathcal{N}(0,\boldsymbol{\Sigma}) \] \[ \Sigma \approx \boldsymbol{U} \boldsymbol{S} \boldsymbol{U}^t = (\boldsymbol{U} \boldsymbol{S}^{1/2} \boldsymbol{U}^t)(\boldsymbol{U} \boldsymbol{S}^{1/2} \boldsymbol{U}^t)^t \]

then

\[ y_{\text{pred}} = (\boldsymbol{U}\, \boldsymbol{S}^{1/2}\,\boldsymbol{U}^t) \times \boldsymbol{Z} \text{ where } Z_i \sim \mathcal{N}(0,1) \]

because \(\boldsymbol{U}^t \, \boldsymbol{U} = I\) since \(\boldsymbol{U}\) is an orthogonal matrix.

\[ n=1000, \quad p=10000 \]

Vecchia’s approximation

Choosing neighbors

Neighbors in practice

Benefits?

• The NN approximation is a huge computational savings over the full likelihood.

  • Dealing with \(n\) \(k \times k\) matrices is much easier than dealing with a single \(n \times n\) matrix.

  • The NN approximation is \(O(nk^3)\) while the full likelihood is \(O(n^3)\).

• The NN approximation is also a huge memory savings over the full likelihood.

• The NNGP is a well-defined spatial process

NNGP Model Results

Performance